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We study the merger of binary neutron stars using different realistic, microphysical nuclear equa¬ 
tions of state, as well as incorporating magnetic field and neutrino cooling effects. In particular, we 
concentrate on the influence of the equation of state on the gravitational wave signature and also 
on its role, in combination with cooling and electromagnetic effects, in determining the properties 
of the hypermassive neutron star resulting from the merger, the production of neutrinos, and the 
characteristics of ejecta from the system. The ejecta we find are consistent with other recent stud¬ 
ies that find soft equations of state produce more ejecta than stiffer equations of state. Moreover, 
the degree of neutron richness increases for softer equations of state. In light of reported kilonova 
observations (associated to GRB 130603B and GRB 060614) and the discovery of relatively low 
abundances of heavy, radioactive elements in deep sea deposits (with respect to possible production 
via supernovae), we speculate that a soft EoS might be preferred—because of its significant produc¬ 
tion of sufficiently neutron rich ejecta—if such events are driven by binary neutron star mergers. 
We also find that realistic magnetic field strengths, obtained with a sub-grid model tuned to capture 
magnetic amplification via the Kelvin-Helmholtz instability at merger, are generally too weak to 
affect the gravitational wave signature post-merger within a time scale of « 10 ms but can have 
subtle effects on the post-merger dynamics. 


Contents 


I. INTRODUCTION 


I. Introduction 0] 

II. Numerical Implementation [2] 

A. Evolution equations [2] 

B. Equations of State and neutrino transport [3] 

C. Analysis [3] 

III. Numerical results 0] 

A. Gravitational Waves [5] 

B. Matter dynamics: Outflow and Ejecta [8] 

C. Neutrino Emission EH 

D. Magnetic Effects [13] 

IV. Conclusions [15] 

A. The primitive solver [17] 

B. Convergence tests [18] 

Acknowledgments 08] 

References [TUI 


Binary neutron star systems (BNS) represent perhaps 
the most exciting and anticipated source for upcoming 
gravitational wave (GW) observations. In contrast to 
other expected binary sources such as neutron star-black 
hole binaries, BNS have been observed and these obser¬ 
vations lead to reasonably robust lower bounds on the 
expected GW detections [I]. In addition, neutron stars, 
being made of matter and supporting strong magnetic 
fields, are more complicated than black holes, especially 
during merger. The coalescence of neutron stars, being a 
highly dynamical violent event, disrupts the stars, pow¬ 
ers electromagnetic events, and produces large neutrino 
bursts. Other effects may produce precursor signals prior 
to merger through either the interaction of stellar mag¬ 
netospheres PHI or crust cracking [5]. 

The merger produces a differentially rotating, hot mas¬ 
sive neutron star (MNS) with a strong neutrino luminos¬ 
ity. Such mergers are suspected engines of short, gamma 
ray bursts (sGRB), arguably the most spectacular elec¬ 
tromagnetic counterpart to a strong GW event. Bina¬ 
ries composed of a black hole and neutron star may also 
power sGRB, but, in order for such binaries to do so, ei¬ 
ther low mass ratios or black holes with significant spin 
are required, e.g. [SHS]. In both systems, the merger 
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may produce a significant neutron rich ejecta which can 
trigger kilonova events m (see also mm)- Remark¬ 
ably, a recent observation of a sGRB led to a follow¬ 
up observation in the infrared indicative of a kilonova 
event HUES]. This observation provides further support 
for the connection between non-vacuum binary mergers 
with sGRBs [161 . 

Moving beyond the broad expectations of BNS mergers 
sketched above, the detailed dynamics of such mergers 
depend on a number of parameters, such as the stellar 
masses, equation of state, and magnetization, to name 
but a few. Fortunately, the combined efforts by the astro¬ 
physics and numerical relativity communities are gradu¬ 
ally refining our understanding of these systems and the 
observational opportunities they present. For instance, 
simulations have: indicated that the merger will give 
rise to a black hole promptly when the total mass of 
the system M tota i > 2 . 8 M© (e.g. [13 EH]); provided an 
approximate range for the final spin of the black hole 
for cases that do collapse 01 ; indicated that the colli¬ 
sion can produce magnetic field strengths in the resulting 
MNS that range from magnetar levels to higher [201122] : 
and contributed gravitational wave templates to aid in 
data analysis efforts (e.g. [23ti25] h 

Most recently, a number of research avenues are at¬ 
tracting considerable interest. A significant effort is ex¬ 
amining possible electromagnetic counterparts, of which 
the aforementioned sGRB are just one example (see 
e.g. [3 [26H30]. Another is the study of binaries with 
more realistic equations of state [23, I3TH33] . In particu¬ 
lar, the hope is that information about the true equation 
of state for NS interiors can be obtained from observa¬ 
tions. Related to such efforts is an intense interest in 
ejecta and outflows [33j|34|, because of their connection 
to possible kilonova events and their role in explaining 
the origin of r-process elements [35]. Neutron star binary 
characteristics, including ejecta, have also been actively 
studied within Newtonian models [3d , where the dynam¬ 
ics sometimes exhibit subtle, but important, differences 
with respect to analogous systems in full general rela¬ 
tivity. Additionally, the role of a magnetic field in BNS 
mergers is also under study. Potentially, an organized 
post-merger field may play a crucial role in the genera¬ 
tion of a jet were a binary to power a sGRB. Since it is 
generally expected that neutron stars have some magne¬ 
tization, determining what electromagnetic effects could 
arise is important. 

Here we explore the dynamics of binary neutron star 
systems with increased realism: we study stars described 
by different realistic equations of state; we account for 
neutrino cooling effects through a leakage algorithm; and 
we consider the role of magnetization in the post-merger 
dynamics. Details of our numerical implementation can 
be found in previous work m eh ezhsu- We summarize 
details in section II, present results in section III, and dis¬ 
cuss these results in section IV. We defer to appendices a 
description of the primitive solver and convergence tests. 


II. NUMERICAL IMPLEMENTATION 

Below we briefly summarize: (i) the evolution equa¬ 
tions for the spacetime and the magnetized fluid; (ii) the 
prescriptions employed for the microphysical equations of 
state and the neutrino cooling scheme; and (iii) a descrip¬ 
tion of the quantities employed to analyze the dynamics. 
Full details of our implementation are described in m- 
We adopt geometrized units where G — c — M 0 = 1 , al¬ 
though some results are reported more naturally in phys¬ 
ical cgs units. 


A. Evolution equations 

The Einstein equations in the presence of both matter 
and radiation can be written as 

G a b = 87r(T a 5 + 7 Z a b), (1) 

where G a b is the Einstein tensor, 7 Z a b is the contribu¬ 
tion from the radiation field, and T a b is the stress energy 
tensor of a magnetized, perfect fluid. 

We solve the Einstein equations by adopting a 3+1 
decomposition in terms of a spacelike foliation. The hy¬ 
persurfaces that constitute this foliation are labeled by a 
time coordinate t with unit normal n a and endowed with 
spatial coordinates x l . We write the spacetime metric as 

ds 2 = —a 2 dt 2 + 7 ij (dx l + ft dt) ( dx J + ft dt ) , (2) 

with a the lapse function, ft the shift vector and 7 ij the 
metric on spatial hypersurfaces. We write the Einstein 
equations in terms of the BSSN-NOK formalism m- 
143] . supplemented with appropriate gauge conditions; the 
“ 1 +log” slicing and the T-driver shift conditions, as de¬ 
tailed in m- In the problem of interest, radiation en¬ 
ergies and stresses are much smaller than matter ener¬ 
gies, or \lZ ab \ < |T a6 |. This observation, together with 
the fact that a neutrino leakage scheme can not possibly 
treat the radiation stress tensor completely consistently, 
allows us to ignore 7 Z a b as a source term for the Einstein 
equations. 

We model the matter as a magnetized perfect fluid 
with a stress energy tensor T a b given by 

Tab = hu a u b + Pg ab + F ac F b ~ \9abF c d,F cd , (3) 

where h is the fluid’s total enthalpy h = p(l + e) + P, and 
{p, e, y e , u a , P} are the rest mass energy density, specific 
internal energy, electron fraction (describing the relative 
abundance of electrons compared to the total number 
of baryons), four-velocity, and pressure of the fluid, re¬ 
spectively. Magnetic effects are included in the Faraday 
tensor, F a b , where a factor l/y/far has been absorbed in 
its definition. Given an equation of state P = P(p, e, Y e ) 
and a relativistic Ohm’s law, the equations determining 
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the magnetized matter dynamics are obtained from con¬ 
servation laws. We adopt the ideal magnetohydrodynam¬ 
ics (MHD) approximation, which states that the fluid is 
described by an isotropic Ohm’s law with perfect con¬ 
ductivity so that the electric field vanishes in the fluid’s 
frame F ab u b = 0, because it provides a simple, yet re¬ 
alistic, model that keeps the number of fields needed to 
describe the electromagnetic effects to a minimum (in 
contrast to, say, a resistive MHD approach [44ll45] h 
The equations of motion consist of the following con¬ 
servation laws 


v a T 6 a = g b , 

V a (T a f,n 6 ) = 0, 

W a (Y e pu a ) = pRy , 

V a *F ab = 0 . 

The sources Q a (= — V c 7^^) and Ry are the radiation 
four-force density and lepton sources which are deter¬ 
mined within the leakage scheme. These equations are 
conservation laws for the stress-energy tensor, matter, 
lepton number, and Maxwell tensor, respectively. Notice 
that, in the absence of lepton source terms (.Ry = 0 ), 
Eq. (| 6 | becomes a conservation law for leptons, similar 
to the baryon conservation law, i.e., Y e is a mass scalar. 
These conservation equations combined with the Einstein 
equations Eq. (!]) comprise the complete set of equations 
we employ to describe the systems of interest. 

The matter equations of motion in Eqs. 0-0 are 
written in balance law form 

d t u + diP (u) = s(u) ( 8 ) 

by defining a set of conservative variables u. The primi¬ 
tive (or physical) quantities are recovered from the con¬ 
servative ones by solving a system of non-linear equa¬ 
tions, described in detail in Appendix |§ A| Due to the 
inability of current fluid codes to evolve very rarefied 
matter, we impose a floor on the density, so that a ten¬ 
uous atmosphere fills the domain outside the stars with 
constant density of 6 x 10 5 g/cm 3 . 

We use finite difference techniques on a regular, Carte¬ 
sian grid to discretize the system das?]. The geomet¬ 
ric fields are discretized using a fourth order accurate 
scheme satisfying the summation by parts rule, while a 
High-Resolution Shock-Capturing method based on the 
HLLE flux formulae with PPM reconstruction are used 
to discretize the fluid and the electromagnetic variables. 
The fluid equations are written in finite difference form 
using a third-order CENO method. The time evolution 
of the resulting equations is performed by using a third 
order accurate Runge-Kutta scheme mam]. To ensure 
sufficient resolution is available in an efficient manner, we 
employ adaptive mesh refinement (AMR) via the HAD 
computational infrastructure that provides distributed, 
Berger-Oliger style AMR |49ll50] with full sub-cycling in 
time, together with an improved treatment of artificial 
boundaries m- 


( 4 ) 

( 5 ) 

(6) 
( 7 ) 


B. Equations of State and neutrino transport 

We use finite-temperature equations of state based 
on several relativistic mean field models of the nu¬ 
clear interaction to model both neutron star matter and 
the hot, lower density material present during and af¬ 
ter the merger. These different nuclear interactions 
give different cold neutron star structures and corre¬ 
spondingly different neutron radii and maximum masses. 
The use of a microphysical EoS over a simple poly¬ 
trope is also needed to model neutrino interactions as 
the neutrino emission, absorption, and scattering rates 
depend sensitively on the matter density, temperature 
and composition. We use publicly available EoS tables 
from www.stellarcollapse.org described in [521. We have 
rewritten some of the library routines for searching the 
table to make them faster and more robust. 

The leakage scheme seeks to account for (1) changes to 
the (electron) lepton number and ( 2 ) the loss of energy 
from the emission of neutrinos. Because the dynamical 
timescale for the post-merger is relatively short, radia¬ 
tion momentum transport and diffusion are expected to 
be subleading effects. As is standard, we consider three 
species of neutrinos: v e for electron neutrinos, v e for elec¬ 
tron antineutrinos, and v x for both tau and muon neu¬ 
trinos and their respective antineutrinos. Our scheme is 
based on the open-source neutrino leakage scheme from 
Ref. [52], also available at www.stellarcollapse.org. In 
particular, the leakage scheme provides the fluid rest 
frame energy sink Q and lepton sink/source Ry due to 
neutrino processes. Ry is the source term for a scalar 
quantity and therefore is the same in all frames. We ex¬ 
press the source term for the energy and momentum in 
an arbitrary frame as 


Ga = Qu a ■ (9) 

Since the effect of neutrino pressure is small [53 and 
difficult to accurately capture with a neutrino leakage 
scheme, we ignore its contribution in the fluid rest frame. 
As mentioned, the details of our implementation can be 
found in Ref. [31], which also describes a novel and effi¬ 
cient calculation of the optical depth. 


C. Analysis 

We compute several quantities to analyze the quality 
of the solution as well as the dynamics and main char¬ 
acteristics of the system. A useful quantity to assess the 
correctness of the numerical solution is the total baryonic 
mass in the domain 


M b = J pW^dx 3 , ( 10 ) 

where W is the Lorentz factor. Because this mass should 
remain constant, we monitor it and also study its con¬ 
vergence in Appendix | § B| 
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Among the three possible signals produced by a 
BNS, the cleanest signal is provided by gravitational 
waves, which are potentially observable by the upcom¬ 
ing aLIGO/VIRGO to a distance of ~ 450 Mpc. We ex¬ 
tract this signal from our simulations by computing the 
Newman-Penrose scalar T 4 , which is directly related to 
the metric perturbation components by 4/4 = — ih x 

in the appropriate Bondi frame [54l [55] (c.f. [56]). We 
compute 4>4 on several spheres centered on the origin and 
with radii of {225, 300, 375}km, i.e, far from the sources, 
in the wave zone. Performing the calculation on multiple 
spheres allows various consistency checks to confirm the 
validity of the signal. As is customary, we expand 4 /4 
in terms of spin-weighted spherical harmonics (with spin 
weight s = — 2 ) 

r* 4 (t, r, 0,<t>)=^2 r) ~ 2 Y ltm (d, 0). (11) 

l,m 


We also extract the instantaneous GW angular frequency 
from the primary (l = m = 2 ) mode through the relation 







2,2 


( 12 ) 


The high cost of BNS simulations makes computing 
the full waveform relevant for matched-filtering analysis 
for data obtained by detectors such as aLIGO/VIRGO 
in full relativity impractical. However, for larger separa¬ 
tions (low frequencies), post-Newtonian (PN) analyses 
provide excellent approximations. Therefore, our sim¬ 
ulations concentrate primarily on the merger stage, 
and hybrid waveforms are constructed by matching the 
simulation waveforms to PN waveforms. Among the 
available PN waveforms, the Taylor-T4 approximation 
has proven quite convenient for binary black hole 
systems E3, and has recently been expanded to account 
for tidal effects [58] (cf. 22)* Our simulations span 
the last 4-6 orbits (depending on the EoS) and can be 
matched to the PN waveforms (by suitably blending 
over a cycle) to construct a complete waveform within 
aLIGO’s frequency window. This allows us to calculate 
the gravitational wave power as a function of frequency 
and contrast it with the expected noise curve of aLIGO. 


Another possible signal is provided by neutrinos 
produced by the shock-heated material during merger 
and afterward. Provided the event takes place rela¬ 
tively nearby, such neutrinos would be observable in 
current detectors such as Super-Kamiokande m and 
IceCube m as well as future detectors such as Hyper- 
Kamiokande [62] . We extract the neutrino luminosity 
and compute from it the distance at which a detection 
could take place. 


Electromagnetic signals are also possible and would 
provide complementary information about the system 
(and with the possibility of covering a much wider set 


of frequencies due to the mature stage of detectors in 
the electromagnetic band). To extract such information, 
reasonable models connecting possible observations with 
particular characteristics of the systems are required. 
The details of such models generally rely on different as¬ 
sumptions and, in many cases, involve longer timescales 
than those that can be reasonably probed within a full 
general relativity simulation. However, one can explore 
the level to which consistent conditions are provided for 
such schemes. In particular, the lifetime of the MNS, 
the timescale for possible black hole formation, and the 
characteristics of any material outflow are key physical 
parameters (see e.g. and references cited therein [29l l63l 
[65]). Thus, we monitor the behavior of the formed 
MNS—e.g. central density and multipolar structure—as 
well as the properties of matter within our simulations, 
supplemented with reasonable approximations to connect 
with possible phenomena at longer timescales. 

We assess the dynamics of all stellar material, at some 
representative times after the merger, by assuming it 
moves along geodesics. We compute its specific energy 
(as measured at infinity) (E = —u t — 1) and consider 
the material unbound if E > 0. We then produce his¬ 
tograms displaying various quantities which characterize 
the state of the bound and unbound material for the dif¬ 
ferent EoS. In particular, to generate a histogram charac¬ 
terizing some quantity x (e.g. temperature), we compute 
a series of histogram heights, hi , for bins extending from 
bi bi +1 (e.g. 10 MeV to 20 MeV) as an integral over 
the whole space hi = f 0(x — &i) 0 (&i+i — x)pdV , where 
0(x) is the Heaviside function. In practice a single sweep 
over the domain suffices because, at each gridpoint, the 
integrand gets added to just one bin. 


III. NUMERICAL RESULTS 

We concentrate on equal-mass binaries with individual 
gravitational masses given by 1.35M©, a value consistent 
with current observations [ 6 d . We consider three distinct 
microphysical EoS: SFHo 67 , DD 2 [53 and NL3 [b 8 j. 
These three EoS produce cold neutron stars whose radii 
cover a large range of values. In particular, the SFHo EoS 
gives the smallest radii (~ 12 km for a mass of 1.35M©) 
among the three and is therefore considered a soft EoS. 
The NL3 EoS yields a large radius (~ 15 km for the 
same mass) and is considered stiff. The DD 2 EoS leads 
to an intermediate radius (« 13 km). We display the 
gravitational mass of isolated neutron stars for each EoS 
in Fig. [l] A horizontal, green line at 1.35M 0 shows the 
stars considered here and its intersection with the three 
curves indicates the three different radii. All three EoS 
can produce neutron stars with a mass of at least 2 M 0 , 
consistent with neutron star mass observations [69] [70]. 

We concentrate on the late stages of the coalescence, 
roughly the last 4-6 orbits (depending on EoS), of equal 
mass, irrotational binary neutron stars. The physical pa¬ 
rameters of the binaries and our grid setup are summa- 
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FIG. 1: Gravitational mass as a function of the circumferen¬ 
tial radius for the three different microphysical EoS considered 
here. The horizontal, light blue line shows the maximum mass 
astrophysically measured for a neutron star Eama and the 
horizontal, green line indicates the mass considered in this pa¬ 
per. These curves were generated with the Magstar solver, 
part of the Lorene package EH- 


rized in Table [TJ To monitor convergence of the solutions, 
we evolved the binary described by the DD2 EoS with 
three different reso lutions obtaining convergent results 
(see Appendix §B). For the two other EoS, we evolved 
only two different resolutions confirming that the differ¬ 
ences between the two were consistent with those of the 
DD2 case. (For reference, unless otherwise specified, we 
set t = 0 as the time when all binaries are separated by 
45km.) 

As mentioned, we include the effects of neutrino cool¬ 
ing and magnetic field in the dynamics. To reduce the 
significant computational overhead that these introduce, 
neutrino cooling is not calculated until a few millisec¬ 
onds before merger when matter heats up and neutrino 
effects become significant. Likewise, we give the stars a 
poloidal magnetic field (with maximum strength 10 13 G) 
a few milliseconds prior to merger. This strategy allows 
us to consider both the role of neutrino production and 
stellar magnetization in an efficient manner because nei¬ 
ther cooling (the stars are cold prior to merger) nor mag¬ 
netic field effects [381 [72] affect the dynamics prior to 
merger. 


A. Gravitational Waves 

At large separations, all binaries follow the same tra¬ 
jectory since the internal structure of the binary con¬ 
stituents only contributes at 5th PN order (see e.g. ES]). 
However, as the binary tightens, differences from a point- 
particle approximation arise and increase as the merger 
ensues. Such differences are intrinsically related to the 
tidal deformability of the stars which is, in turn, directly 


tied to the EoS. 

We summarize the pre-merger dynamics of our three 
binaries in Fig. |2j The (coordinate) separation, the GW 
signal \Eg (l = 2,m = 2 component), and the instanta¬ 
neous GW frequency, all as functions of time, are shown. 
Also included in the bottom frame is the gravitational 
wave frequency as predicted by the T4 PN approxima¬ 
tion [57]. During the first eight milliseconds, all the bi¬ 
naries agree with each other and with the T4 approxi¬ 
mation. Only at times closer to merger do the binaries 
begin to differentiate from each other and from the point- 
particle approximation. 

The first EoS to show differences is NL3, followed by 
DD2, and then SFHo. This is the expected behavior 
as the NL3 EoS corresponds to the largest neutron star 
radii, DD2 the intermediate and SFHo the smallest radii. 
Tidal effects are proportional to fe-R 5 (where &2 is the 
tidal Love number and R is the stellar radius), and there¬ 
fore, for fixed mass, differences in radii are dominant. 
We take advantage of this early agreement to compute a 
complete wavetrain by creating a hybrid of the PN and 
our simulated waveforms which we show in Fig. [6] and 
discuss later. Of course, the merger itself is highly non¬ 
linear, and so an examination of the merger waveforms 
does not involve any further comparison with PN results. 
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FIG. 2: Orbital parameters displayed as a function of time 
for the three different EoS considered here. (Top) Separation 
between the stars, computed as the coordinate distance be¬ 
tween the maximums of the density. (Middle) Primary grav¬ 
itational wave mode l — m — 2. (Bottom) The freque ncy of 
the primary GW mode where /gw = w / (27r) (see Eq. (12)). 
Also shown is the frequency obtained from the Taylor T4 
approximation. Note that differences among the EoS only 
appear at times near merger. 


Considering now the merger, we shift times in each case 
such that t = 0 refers to the moment when the maximum 
density first occurs at the origin of the binary. Fig. [3] 
presents the gravitational waveforms with time measured 
from the merger. Note that all three signals prior to 
merger (negative times) oscillate within a common am- 
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EoS 

M o adm [M 0 ] 

Jo DM [GMq/c] 

nA 1 ) [Mq] 

[km] 

D 0 [rad/s] 

/o GW [Hz] 

Afeject [A/q] 

Ax m i n [m] 

NL3 

2.7 

7.40 

1.470 

12.63 

1778 

566 

1.6 x 10 “ e 

230 

DD2 

2.7 

7.39 

1.485 

10.88 

1776 

565 

4.3 x 10 -4 

230 (275 & 192) 

SFHo 

2.7 

7.38 

1.498 

9.47 

1775 

565 

3.2 x 10 _a 

230 


TABLE I: Summary of the parameters of the binary systems considered in this work. The initial data were computed using 
the Bin star solver from the Lorene package m by imposing beta-equilibrium and a constant temperature T = 0.02 MeV. 
All the binaries start from an initial separation of 45 km and the outer boundary is located at 750 km: the ADM mass M^ dm 
and angular momentum J^ DM of the system, the baryon mass of each star and its coordinate radius r^ l \ the initial orbital 
angular frequency Do of the system, the initial GW frequency = Qo/tt and the best resolution Ax m i n coveri ng b oth the 
stars. For the DD2 EoS, higher and lower resolutions were carried out for the study of convergence (see Appendix §B). 


plitude envelope. This behavior can be understood with 
an effective, and approximate, one-body problem point 
of view. In such a problem, the central object has a mass 
approximately given by the total mass of the system and 
radiative properties can be captured by the infall of a 
particle of reduced mass fi = m\/(mi 4 - 777 , 2 ). Shifting 
time so that all cases merge at t = 0 allows for identi¬ 
fying a natural common time, and phase differences are 
a consequence of the particle being released at different 
times with respect to the shifted time. Such an approach 
has been at the core of a number of effective prescriptions 
for capturing various aspects of the two body problem in 
general relativity (e.g. [74HZZ]). 

Fig. [3] also allows for a clear comparison of the strength 
of the gravitational waves and the waveform periods. The 
SFHo EoS achieves the largest amplitude during merger 
as a natural consequence of yielding the smallest stel¬ 
lar radii. A smaller radius means that the merger takes 
place deeper in the gravitational potential of the binary, 
and hence more energy is released, resulting in a more 
dynamical and violent collision. In contrast, the weakest 
peak-signal with longest period is produced by the NL3 
EoS with its large radii while the DD2 binary is interme¬ 
diate between the two extremes considered. 

Qualitatively, the mergers for all three EoS proceed 
similarly. The stellar orbit tightens as GW are emitted, 
and, as the stars get closer, tidal deformations grow. Ev¬ 
idence for such deformations can be appreciated, for the 
DD2 case, in the first two frames of Fig. [ 4 ] (We cau¬ 
tion however that coordinate effects are not taken into 
account in the figure; observable effects of tidal defor¬ 
mations are better represented in the gravitational wave¬ 
forms obtained.) After the merger, strong radial oscilla¬ 
tions in the remnant lead to peak densities bouncing out¬ 
ward (with respect to the co-rotating frame), as shown 
in the last four frames of Fig. [4] These bounces are also 
apparent in the separation plot in Fig. |2j The oscil¬ 
lations however are stronger—larger in amplitude and 
longer lasting—as the EoS becomes stiffer, with mergers 
at earlier times and smaller angular velocities (as can be 
appreciated in Figs. [ 2 ] and [3| . 

The resulting hot MNS will cool down as well as lose 
energy and angular momentum via emission of neutrinos 
and gravitational waves. Furthermore, its velocity dis¬ 
tribution will be redistributed through angular momen¬ 
tum transport by hydrodynamical and electromagnetic 


effects. Collapse to a black hole will eventually take place 
if the mass of the MNS is higher than the maximum al¬ 
lowed mass for a given EoS. Inspection of Fig. [TJ together 
with the observation that the amount of unbound mass 
is relatively small in all cases, indicates that at least the 
SFHo case is expected to collapse to a black hole at some 
point (see also [T8l [33] ). Such an expectation arises nat¬ 
urally because the maximum mass allowed with this EoS 
for a non-rotating star is ~ 2.05M©, and, while rotation 
helps increase the maximum allowable mass, it can not 
do so beyond « 20% m- Indeed, for the times studied 
in this work, the SFHo binary collapses to a black hole 
~ 7 ms after merger. 



FIG. 3: Comparison of the near-merger behavior of the GW 
signal, as indicated by the primary mode of T 4 . The different 
signals have been shifted in time such that t — 0 corresponds 
to the time when the separation between the density maxi- 
mums vanishes. The significant differences among the signals 
near merger indicate that information about the internal stel¬ 
lar structure (e.g. details of the EoS, radius, compactness) is 
encoded within the gravitational wave signature. 

Additional information about the properties of the 
neutron star are provided by the GW after the merger. 
Fig. [5] displays the power spectral density of the Fourier 
transform for the main GW mode at this stage. The 
spectrum shows several characteristic peaks. The dom¬ 
inant peak (the one with the most power), / pea k, is as¬ 
sociated with the rotational behavior and quadrupolar 
structure of the MNS. As mentioned, a softer EoS re- 
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t= 13.1 ms t = 13.7 ms 




FIG. 4: Density on the equatorial plane for the DD2 EoS. 
The stars merge at t ~ 14.3 ms. The first two frames show the 
system pre-merger as tidal deformations arise. Also visible in 
the final four frames is the rebound of the stellar material 
post-merger. 

suits in a merger at higher frequency as the smaller radius 
stars coalesce. This effect is apparent when comparing 
the spectra for the different EoS in Fig. [5| The peak 
frequencies are also presented in Table [II| 

Beyond this qualitative discussion, the peak frequency 
is slightly higher than twice the orbital frequency at the 
time of merger. Since there will be further compression of 
the merged MNS, its rotational frequency will be higher 
than the frequency estimated by Kepler’s law at the mo¬ 
ment of contact. Nevertheless, the Keplerian orbital fre¬ 
quency can be used to estimate / pea k as jl 12] 

/peak//peak « {Rj/Ri? /2 , (13) 

where i and j denote different binaries with stellar radii 
given by Ri and Rj , respectively. Interestingly, several 
other subdominant peaks exist at lower frequencies that 
arise from oscillations of the MNS. All these peaks offer 
interesting future detectability prospects for 3rd genera¬ 
tion detectors (aLIGO/VIRGO/KAGRA) that will have 
specifically tuned configurations to increase sensitivity in 
the higher frequency range or even signal-stacking possi¬ 
bilities (e.g. US])- ' 

Recently, Ref. [80] indicated that these low frequency, 
secondary modes can be traced to two different mech¬ 
anisms, although this is still a matter of debate (see 


EoS 

h 

h 

h 

u 

M/R 

/peak 

/spiral 

0 

1 

NL3 

2.0 

1.5 

0.8 

- 

0.135 

2.2 

1.6 

1.2 

DD2 

2.3 

2.0 

1.8 

1.62 

0.15 

2.6 

1.9 

1.5 

SFHo 

3.4 

2.6 

2.2 

1.6 

0.167 

3.2 

2.4 

2.1 


TABLE II: Post-merger oscillation frequencies in kHz. The 
frequencies of the various peaks of the post-merger GW spec¬ 
trum (see Fig. [ 5 ]) are shown as /i,/ 2 ,/ 3 , and The com¬ 
pactness of the neutron stars for each EoS (obtained from the 
initial data for an isolated star) along with the predicted peak 
frequencies from Ref [80] . The decent agreement between /1 
and /peak, between / 2 and /spiral, and between either /3 or 
/4 with /2-0 suggests consistency with the model presented 
in [80]. 

for instance [ED EH). Using the same naming conven¬ 
tions of [80] . the first, secondary peak / 2-0 corresponds 
to a non-linear combination of the dominant quadrupo- 
lar mode with the quasi-radial oscillations of the star. 
The other sub-dominant peak, / spira i, is produced by 
a strong spiral deformation induced at the merger time 
that lasts for just a few rotational periods. We calculate 
these frequencies for our simulations using the estimates 
in Ref. [80] . which are based on the compactness of the 
neutron stars, and include these in Table [Tlj Comparing 
the frequencies extracted from our simulations with the 
estimates, we find that the simulation values for both 
the dominant, / 1 , and secondary, / 2 , peaks are consis¬ 
tent with / P eak and / sp i ra L respectively. Also, the sim¬ 
ulation frequencies {/ 3 ,/ 4 } are in reasonable agreement 
with /2-0 for the softest and intermediate EoS (SFHo 
and DD2). 



FIG. 5: Power spectral densities (PSD) of the gravitational 
waveforms after the merger. The Fourier transforms display 
the main oscillation modes of the MNS. The arrows indicate 
particular peaks with their corresponding frequencies R pre¬ 
sented in Table mi 

Finally, we return to the complete gravitational wave- 
train. As discussed, the early behavior is well approxi¬ 
mated by the PN (Taylor T4) expansion which can be 






























augmented to include leading order tidal effects [58H59]. 
We thus obtain a more complete waveform by matching 
this augmented PN signal at early times to our numeri¬ 
cal signal before merger to form a hybrid waveform. The 
matching is accomplished by choosing a cycle after initial 
numerical transients but still well before merger. Within 
this cycle, a weighted average of the two waveforms is 
calculated using a tanh function that transitions from 
0 to 1. We show (a portion of) the hybrid waveforms 
resulting from the three EoS in Fig. [6] 

Additionally, we compare the strength of the signals 
to the expected noise curve of Advanced LIGO as illus¬ 
trated in Fig. [7| which presents the Fourier spectrum of 
the GW (multiplied by y/f) for binaries at 100 Mpc ver¬ 
sus frequency together with two fits to the aLIGO Noise 
curve (the so called “No Signal Recycling mirror” and the 
“zero-detuned high-power case”). As mentioned, differ¬ 
ences among the waveforms arise only at high frequencies 
when tidal effects become relevant. The NL3 waveform 
departs first from the other two cases because the larger 
stellar radii lead to stronger tidal effects. The differences 
among the three EoS arise below the aLIGO noise-curve 
in the case without a signal recycling mirror but could be 
captured in the “broadband” (zero-detuned, high power 
case) configuration. 

Finally, we can compute the “distinguishability” be¬ 
tween the different waveforms (as described in e.g. [32] . 
see also |Sil2J[S5) using the “zero-detuned high-power 
case” noise curve in aLIGO. Without allowing for fre¬ 
quency shifts (and not normalizing our waveforms), we 
find the values for Sh = (hi — hj , hi — hj) are given by 
~ {1.33, 1.36, 1.45} when considering the pairs SFHo- 
DD2, SFHo-NL3 and DD2-NL3 respectively within the 
frequency window / E [300, 8000] Hz. Since all cases sat¬ 
isfy Sh > 1 the gravitational waves will be distinguish¬ 
able from each other by aLIGO with the “zero-detuned 
high-power case” configuration. In contrast, these val¬ 
ues reduce to Sh ~ {0.27, 0.28, 0.32} for the “No Signal 
Recycling mirror” configuration. 


B. Matter dynamics: Outflow and Ejecta 

Having presented characteristics of the GW produced 
mainly by the bulk motion of the stars, we now turn to 
an analysis of the material outflow resulting from the col¬ 
lision. While such motion leads only to subleading effects 
on the bulk dynamics (and hence on GW emission), the 
details of any outflow can be very important for produc¬ 
ing electromagnetic and neutrino signals. 

In particular, among the outflow, material bound to 
the system either remains part of the resulting MNS or 
otherwise contributes to an accretion disk which may po¬ 
tentially power a sGRB (see e.g. [86]). On long time 
scales, such a disk can induce winds that produce op¬ 
tical electromagnetic emission via the radioactive decay 
of r-process elements [13]. Furthermore, the fallback of 
bound material may trigger the eventual collapse of the 



FIG. 6: Hybrid waveforms obtained by matching the PN 
waveform with the respective EoS waveform from the evo¬ 
lution of the binary. These waveforms are used for comparing 
signals with aLIGO noise curves in Fig. [7] 



f (Hz) 


FIG. 7: Fourier spectrum of the GW signals for the three bina¬ 
ries. The dotted and dashed-dotted monotonic curves at high 
frequencies show two fits to the noise power spectrum \fSh of 
aLIGO (specifically the “zero-detuned high-power” and the 
“No Signal Recycling Mirror” (NSRM) cases, see [85]). 


remnant star to a black hole which, in turn, could power 
an intense burst of electromagnetic radiation [26] [87]. 

Another particularly interesting mechanism powered 
by outflows is the proposed kilonova in which material 
ejected by the collision undergoes r-process nucleosynthe¬ 
sis. Being ejected from the core region, the behavior of 
this material and its consequences can be analyzed with¬ 
out worrying about the rather involved behavior of the 
central region. Indeed, as discussed in Pi, a promising 
source of electromagnetic signals is the radioactive de¬ 
cay of r-process elements formed in the ejecta (see also, 
e.g. [12]). Such signals have specific characteristics ob¬ 
servable on the order of days after the merger mm- 
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FIG. 8: Density, temperature, electron fraction and neutrino 
luminosity rate for the different EoS in the z — 0 plane at 
roughly t — 3 ms after the merger. Shown are the NL3 EoS 
(left), DD2 EoS (middle), and SFHo (right). 


FIG. 9: Density, temperature, electron fraction, and neutrino 
luminosity rate for the different EoS in the y — 0 plane at 
roughly t — 3 ms after the merger. Shown are the NL3 EoS 
(left), DD2 EoS (middle), and SFHo (right). 


Importantly, such characteristics can have a rather tight 
dependence on the parameters of the binary and the EoS 
describing the neutron stars involved. 

Given the importance of the outflow, in the following 
we compare the properties of both bound and unbound 
material among the three EoS. 

We begin our analysis of the merger and postmerger 
by presenting snapshots of important quantities just a 
few milliseconds after merger. Figs. [8] and [9] display the 
density, temperature, electron-fraction, and neutrino lu¬ 
minosity along the z = 0 and y = 0 planes for the three 
EoS. The density snapshots show clear differences among 
the EOS, with the SFHo having a much more distributed, 
high density region than the NL3 remnant. This differ¬ 
ence can be understood in terms of the stiffness of the 
EoS. In particular, the SFHo, being the softest of these 
EoS with the most compact stars, results in a more vi¬ 
olent merger; the compact stars have further to travel 
to merge and have higher impact velocities. The vio¬ 
lence serves to distribute the stellar material more than 
is seen with the stiffer EoS and, in particular, yields more 
ejecta. Likewise, the temperature of the material is gen¬ 
erally higher for the softer EoS and, as a consequence, 
the neutrino production is also higher (we discuss fur¬ 
ther neutrino-related aspects in the next subsection). 

To be more quantitative, we present a number of 


histograms that characterize the statistics of this post¬ 
merger material, measured at roughly t — 3 ms after the 
merger. We concentrate on the behavior of the matter at 
this early time to ensure that atmosphere effects do not 
contaminate the extracted values of the outflows. Note 
however, that this implies that our quantities of unbound 
material might be consequently lower than the actual val¬ 
ues, as matter can be ejected not only due to the tidal 
interactions, but also because of thermal pressure, which 
increases significantly at the merger due to shock heating. 
We present the velocity distribution in Fig. [lOj as well as 
directional dependence (as measured with respect to the 


direction of the orbital angular momentum) in Fig. 11 


the temperature in Fig. |T2j and the electron-fraction in 
Fig. [13] An important distinction is whether this mate¬ 
rial is bound to the remnant or whether it is unbound 
ejecta. 

When looking at the unbound material, the reader 
should understand that the histograms display the frac¬ 
tional masses of material with some property. The huge 
disparity in unbound material among the EoS is therefore 
not immediately apparent. Therefore, for the unbound 
matter we include histograms of the absolute quantities 
in the insets; the stiff EoS material is barely visible. 
In particular, SFHo has 0.0032M© unbound material, 
the DD2 has roughly 13% of this amount and the NL3 
has only 0.05% (these amounts are also included in Ta- 
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v/c v/c 


FIG. 10: Distributions of speeds for stellar material. (Left) 
The fraction of ejected mass binned according to velocity. 
The inset shows instead the absolute mass binned in the same 
way. The inset demonstrates that the NL3 case has essentially 
no unbound mass compared to either the DD2 or SFHo cases, 
and the SFHo has significantly more unbound mass than the 
DD2 case. (Right) The fraction of bound material with some 
particular velocity. 



polar angle [radians] 


polar angle [radians] 


FIG. 11: Polar distributions of material for the different EoS 
cases. Stellar material is binned according to its polar angle. 
(Left) Ejected material. (Right) Bound material. Both 
bound and ejected material is mostly found near the orbital 
plane. The asymmetry across the orbital plane for the NL3 
unbound mass is very pronounced and probably indicative 
roughly of our computational error because the amount of 
unbound material for that case is so small. 


gions. Since such winds can also trigger r-processes 
though seemingly with a larger value of Y e mm] and 
thus likely to induce signals through radioactive decay in 
the optical and ultra-violet bands. 


Fig. [T2| contrasts the temperature distribution among 
the EoS. The temperatures of unbound material are 
definitively low for all three EoS. However, the bound 
material of the SFHo case is generally much hotter than 
the other two EoS (i.e., with a distribution peaked at 
temperatures in (10-20) MeV in the SFHo case as op¬ 
posed to 0-5 MeV for the stiffer ones, as well as a tail with 
temperatures higher by lOMeV with respect to the tail 
in temperatures measured for the DD2 and NL3 cases. 
As mentioned, the soft EoS yields a more violent colli¬ 
sion with concomitant higher temperatures. The largest 
fraction of SFHo bound material is close to 30 MeV while 
the stiff peaks at roughly 10-20 MeV. 

Finally, Fig. 13 shows the electron fraction distribution 
and, in particular, that the unbound material is gener¬ 
ally quite neutron rich. The calculated values have a 
significant fraction below Y e < 0.22, a rather robust con¬ 
dition for r-process nucleosynthesis being able to create 
high mass-number elements m j89l . Consequently, an 
electromagnetic signal peaking in the infrared would be 
expected in such cases. 



temperature [MeV] temperature [MeV] 


ble[T|. That a soft EoS yields more ejecta agrees with 
recent results of Ref. [33], and correlates with the colli¬ 
sion occurring deeper in the gravitational potential than 
for the stiffer EoS. The observation of an apparent kilo- 
nova [14, 15], which would require significant ejecta, is 
therefore suggestive of a soft-intermediate EoS, provided 
the binary is composed of two neutron stars with similar 
masses. We will return to this point in the conclusion. 

Fig. [lO] illustrates the fraction of ejected and bound 
mass, binned according to velocity. While the distribu¬ 
tion in speeds of the bound material is within a similar 
range v G (0,0.5) c, although the velocities of the ejecta 
are quite different given the narrower peak for the stiffer 
equation of state. The stiffest EoS, NL3, has large veloc¬ 
ities. It is not clear, however, if these high velocities have 
any physical significance because the quantity of ejecta 
is so tiny. 

We display the polar angular distribution of the post¬ 
merger material in Fig. [TT] The bound material is mostly 
in the equatorial plane, regardless of EoS. Unbound ma¬ 
terial is similarly distributed mostly near the equato¬ 
rial plane. Both these observations imply the system 
avoids baryon poisoning along the polar regions and so 
winds could be preferentially channeled along such re- 


FIG. 12: Distributions in temperature of stellar material. 
(Left) Unbound material is generally cool for all three EoS. 
(Right) Bound material. The temperature of the SFHo case 
is noticeably higher than that of the stiff EoS. 



FIG. 13: Distributions of electron fraction for the stellar ma¬ 
terial. (Left) Unbound material for the different EoS. The 
inset shows the unrescaled amount of unbound material with 
the given electron fraction. Note that the outflow material is 
neutron rich (Y e < 0.3). (Right) Bound material. The in¬ 
set here does not show absolute quantities, but instead shows 
the same information as the main plot but with much smaller 
vertical scales. 
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C. Neutrino Emission 

We compute neutrino cooling and the corresponding 
luminosity through a leakage scheme as reported in m • 
Fig. p| shows the luminosities for each neutrino species 
and the total neutrino luminosity as a function of time 
for each EoS. Following along the trend of the previous 
results which showed the largest gravitational wave emis¬ 
sion and largest ejecta mass for the SFHo EoS, we also 
see the largest neutrino luminosity for this EoS. This re¬ 
flects the fact that the temperature and decompressed 
mass achieved in this case is the highest. Irrespective 
of the particular EoS, it is evident that the luminosi¬ 
ties for both electron neutrinos and electron antineutri¬ 
nos become roughly comparable while the heavy-lepton 
neutrino types are less luminous (we note that the heavy 
lepton luminosity shown here is actually the contribution 
for all four heavy-lepton species). This dominance has al¬ 
ready been observed in other binary neutron star merg¬ 
ers mm i90' and also neutron star-black hole merg¬ 
ers m, and in simulations with more complete neutrino 
transport m E2U- This phenomena is due to neutron 
rich material being shock heated and decompressed [92] . 
This neutron rich, low density, hot material is initially far 
below the new /3-equilibrium value of Y e and will prefer¬ 
entially emit electron antineutrinos until equilibrium is 
reestablished. 



a 


a* 


o* 



2 4 

time [ms] 


FIG. 14: Neutrino luminosities for the different EoS. 


Given the high neutrino luminosities achieved during 
the merger it is interesting to study their possible de¬ 
tection in future and current facilities. For the pur¬ 
poses of detection on Earth, we are most interested in 
the spectrum and luminosity of the electron antineutri¬ 
nos, as these neutrinos will produce the largest signal 
in Earth-based water Cherenkov detectors, such as the 
current Super-Kamiokande or IceCube detectors, or the 
future Hyper-Kamiokande detector. It is true that neu¬ 
trino flavor oscillations can occur due to one or more 


of the following: (1) neutrino-neutrino interactions, (2) 
neutrino-matter interactions, or (3) oscillations in vac¬ 
uum. In these cases, some fraction of the electron an¬ 
tineutrino spectrum will oscillate into heavy lepton an¬ 
tineutrinos and vice versa. The precise influence of neu¬ 
trino oscillations on the electron antineutrino spectrum 
observed in Earth based detectors is unknown, but de¬ 
pends sensitively on the unknown neutrino mass hierar¬ 
chy and the neutrino and matter conditions in the merger 
remnant [93 . For this reason, in our analysis we neglect 
oscillations and assume the electron antineutrino signal 
at Earth is the emitted signal from our merger calcula¬ 
tion with a gravitational redshift. We note that Ref. m 
recently showed that the electron antineutrino luminos¬ 
ity predicted from their leakage scheme, which is similar 
in design to ours, matches the prediction from a more 
complete neutrino transport calculation to better than 
~ 30%. Furthermore, the neutrino luminosities shown 
in Fig. 14 agree well with the recent results of 


who 

study similar EoS and initial conditions. The leakage 
scheme also predicts a lepton loss rate that, in theory, 
provides an average energy for the emitted neutrinos. 
However, this average energy prediction is typically too 
high because it includes contributions from very high en¬ 
ergy neutrinos leaking out from optically thick zones. In 
reality these neutrinos thermalize as they diffuse out. A 
more representative average energy is one determined by 
the thermodynamic conditions near the neutrinosphere, 
the location where the neutrinos leave thermodynamic 
equilibrium with the matter. 


To make a more accurate determination of this neu¬ 
trino average energy, we use the formalism described in 
Ref. [94]. We chose to examine the disk as seen by an 
observer looking down from z = Too, i.e. looking face- 
on towards the disk. To locate the neutrinosphere of the 
remnant disks for this direction, rays are traced down 
onto the disk from z = +oo and halted when the in¬ 
tegrated optical depth reaches 2/3. We repeat this for 
many (x,y) pairs. We ignore general relativistic effects 
on the ray propagation. The stopping location sets the 
neutrinosphere surface, <S. During the integration, at a 
given spatial location, we assume the opacity is the av¬ 
erage over a Fermi-Dirac, zero chemical potential neu¬ 
trino number distribution with a temperature equal to 
the local matter temperature. As contributions to the 
opacity we include: (1) elastic scattering of all neutrino 
species on neutrons, protons, and electrons; (2) neutrino 
capture on protons (neutrons) for electron antineutrinos 
(electron neutrinos); and (3) neutrino-antineutrino anni¬ 
hilation. As we include scattering cross sections, our neu¬ 
trino surfaces correspond to the surface of last scattering, 
not necessarily the surface where the neutrinos thermo¬ 
dynamically decouple. The differences between these two 
surfaces is small for electron-type neutrinos and antineu¬ 
trinos because the cross sections are dominated by the 
charged current processes, but this difference is larger 
for the heavy lepton neutrinos which do not have such 
interactions in BNS merger remnants. That being said, 
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to determine the average energy of neutrinos emitted in 
this direction, for a particular neutrino species, we aver¬ 
age over <S, and compute 


(E vt ) = L dA I dEE<j) Vi (E,S) / j dA J dEfo^EtS), 

(14) 

where the integral is carried out over the neutrinosphere 
surface <S, determined by the rays; the Fermi-Dirac flux, 
ct>(E,S ), is given by 

= < 15 > 


where /fd is the Fermi-Dirac distribution with a tem¬ 
perature equal to the local matter temperature on the 
surface S. The neutrino energy will be gravitationally 
redshifted as neutrinos stream away from the merger rem¬ 
nant. We estimate the neutrino energy at infinity by 
transforming the emitted quantities using the invariance 
of phase-space density, 

1 dN _ / FD I 
c 2 d 3 xd 3 p ~ h 3 c 2 ~ E 3 ’ [ } 


where / is the specific intensity, and 

dN _ c 2 dN 

d 3 xd 3 p ~ WdAdtdEdn' 


(17) 


To proceed, let us denote the quantities measured by an 
observer with a “tilde”, (i.e. the observed energy is E , the 
area differential is cL4, etc.) while the analogous emitted 
quantities are E and dA, etc. Now, by invariance of I/E 3 
we have 


1 dN _ 1 dN 
c 2 d 3 xd 3 p c 2 d 3 xd 3 p 

Focusing on the left hand side, 

1 dN _ 1 dN _ / FD 
c 2 d 3 xd 3 p ~ E 2 dEdAdtdn ~ h 3 c 2 


(18) 


(19) 


we obtain, 


d A= [ E^E 2 dEdAdn. (20) 

dt J h 6 c 2 

We can now rewrite Eq. ( |20| ) in terms of the emitted 
quantities given that E = ^gdoE, assuming that the 
emission is isotropic, and that distances are stretched 
by a factor of (1 — r s /r) -1 / 2 = g (utilizing the 
Schwarzschild metric for simplicity). Thus, Eq. pQ| ) can 
be re-expressed as 


dN f /fd 3/2 r 2 7 T?dA f /fd 1/2^2 1^1 a 

S = J E dE » = b J E dEdA ’ 

( 21 ) 

and thus 

^ = j glNiE,S)dEdA. ( 22 ) 


Time 

(EpA 

{Eu e ) 


R v 

[ms] 

[MeV] 

[MeV] 

[10 53 erg/s] 

[#/ms] 

2.5 (NL3) 

18.5 (22.4) 

15.2 (18.3) 

0.71 

18.1 

3.0 (DD2) 

18.3 (22.1) 

14.6 (17.4) 

1.1 

28.2 

3.2 (SFHo) 

24.6 (29.7) 

23.5 (28.3) 

3.5 

120.8 

8.4 (NL3) 

13.4 (15.6) 

9.8 (11.3) 

1.1 

20.7 

7.9 (DD2) 

13.2 (16.1) 

10.2 (12.4) 

1.6 

29.6 


TABLE III: Observed (emitted, with no gravitational redshift 
taken into account in parenthesis) neutrino average energies 
for electron antineutrinos and electron neutrinos, the electron 
antineutrino luminosity predicted from the leakage scheme 
(from Fig. 14), and the estimated instantaneous detection rate 
in Super-Kamiokande for a merger at lOkpc from Earth via 
Eq. [25] for two different times and three EoS. In the text, we 
give estimates of the integrated event rate up to the 6 ms after 
merger. 


We can proceed in a similar way for EdN/dt = dE/dt , 
and find, 

(1 § = J + \ dEdA = j g 00 cf>(E,S)EdEdA. 

(23) 

This agrees with the well known result that the lumi¬ 
nosity at infinity is = gooL in the case where the 
neutrino surface is at a constant redshift. 

Eq. ( [23] ) contains goo which must be evaluated at 
the neutrino surface. Then the average observed en¬ 
ergy is ( E Ui ) = (dE / dt) / (dN / dt ). For simplicity, when 
evaluating the numerical value of goo, we employ the 
Schwarzschild metric (in Schwarzschild coordinates) with 
a central black hole mass of 2.7 M©. Inside 25 km (^3 
Schwarzschild radii), to avoid spuriously overestimating 
the amount of redshift, we fix the value of the metric to 
the value at 25 km, i.e. |goo|( r < 25 km) ~ 0.68. (For the 
purpose of our estimates, this simple evaluation does not 
depend sensitively on these details however.) We present 
the average energies computed this way for several times 
and EoS in Table IIII1 

In Fig. [l5j we show snapshots of the spatial distribu¬ 
tion of the temperature at the neutrinosphere at ~ 2.5- 
3 ms after merger for the three EoS and for each electron- 
type neutrino species. As the matter is primarily neutron 
rich, the electron neutrino neutrinosphere has a larger 
extent than the electron antineutrinosphere. The SFHo 
merger produces the most decompressed material and 
ejecta and the highest temperatures; thus powering the 
largest neutrino luminosity. This EoS also has larger 
and generally hotter neutrinospheres for all species. For 
the NL3 and DD2 EoS, the neutrinosphere is generally 
confined to the MNS except for the region around the 
tidal tails. These tails are cold and have little mass (see 
§III.B). 

Armed with the neutrino luminosities and average en¬ 
ergies, we can estimate the detection rate at Earth-based 
neutrino detectors. Given the uncertainties introduced 
by the leakage scheme, we make the following assump- 
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FIG. 15: Color map of the temperature at the scattering 
neutrinospheres at ~ 2.5-3 ms after merger for the NL3 EoS 
(left), DD2 EoS (center) and SFHo EoS (right). The top pan¬ 
els show the electron neutrinosphere, and the bottom panels 
show the electron antineutrinosphere. 


tions: (1) We assume that the luminosities shown in 
Fig. [14] are radiated isotropically from the merger. Ob¬ 
servers situated along the merger axis would observe a 
larger than average neutrino luminosity, while observers 
situated on, or near, the equator would observe smaller 
average luminosities [95]. Therefore, the following pre¬ 
dictions are the average expected rate over many sky 
orientations. (2) We assume the spectrum of neutrinos 
follows a zero chemical potential, Fermi-Dirac spectrum 
with a neutrino temperature that gives the average en¬ 
ergy computed above ((E Ui ) ~ 3.157k.). (3) We assume 
that if an interaction occurs in the volume of the de¬ 
tector, it will be registered as an event. This may not 
be true for the lowest energy neutrinos, but such events 
should be few in number due to the low cross section. (4) 
Finally, as mentioned above, we neglect neutrino oscilla¬ 
tions which may alter the electron antineutrino signal at 
Earth. Given these assumptions, the detection rate in a 
water-Cherenkov detector is then roughly given as 


R v ~ N t x n v x 


I J%„<riBDdE 
Jjf TdE 


(24) 


where Nt is the number of target protons, which 
for the 32kT Super-Kamiokande is equal to 2/18 x 
32 x 10 9 g/m amu ~ 2.1 x 10 33 , n„ = L/{E„)/(4nD 2 ) 
is the total number-flux of neutrinos a distance D 
from the source, and f aiBr>dE/ f E% dE is the 
spectrum weighted average cross section with E\ = 
E 2 /[exp(E/T) + 1] equal to the zero-chemical poten¬ 
tial Fermi-Dirac neutrino number distribution and where 
ctibd is the inverse-beta decay cross section, which can 
be approximated as ctibd a 0 /(4m 2 e )(l + 3g\)E% with 
(Jo = 1-761 x 10 -44 cm 2 , m e = 0.511 MeV, g A = -1.254. 


All together, the detection rate is 


Rv 


21.1 


' 32kT ' 

r l v I 

r (Eu) ] 

10 kpc 

_ AT wa ter _ 

10 53 erg/s 

15 MeV _ 

D 


(25) 

We have used Eq. [25] to determine detection rates for 
the times and EoS listed in Table hd We take the elec¬ 
tron antineutrino luminosities from Fig. [14] for the cor¬ 
responding time and EoS (which we include in Table III 
for completeness). The resulting instantaneous detection 
rates for a merger located at 10 kpc from Earth in a detec¬ 
tor like Super-Kamiokande are listed in the last column of 
Table Hill To estimate the total number of neutrinos that 
can be detected, we take a fiducial length of time of 6 ms 
and integrate over Eq. [25] using the electron antineutrino 


luminosities shown in Fig. 14 and assuming a constant 
(Ep e ) equal to the value near 3ms in Table III For the 
NL3, DD2, and SFHo EoS, we calculate ^135, «176, 
~405 events, respectively, in the first 6 ms after merger 
in a Super-Kamiokande like detector and a merger event 
at ^10 kpc. The merger with the SFHo EoS predicts 
much higher rates because both the luminosity and the 
average energy of the neutrinos is larger due to the higher 
matter temperatures achieved in the merger. If the re¬ 
sulting MNS collapses to a black hole (as in the SFHo 
case studied here), this could lead to a sharp reduction 
in the neutrino luminosity (e.g. Ref. [33]), otherwise, we 
would expect a slower decline in the neutrino luminos¬ 
ity as the MNS and disk cool. Ultimately many more 
events could be detected in the MNS/disk cooling phase, 
and both the total number, and signal duration could 
strongly depend on the EoS [96] . 

One can ask how far would we be able to detect a BNS 
merger in neutrinos. As the distance of the merger in¬ 
creases, the rate of detection decreases quadratically, as 
seen in Eq. (25). For mergers similar to the ones studied 


here, but located in the Andromeda galaxy ^800 kpc 
from Earth rather than 10 kpc, the detection rates are 
^6,400 times smaller. Consequently, we would not ex¬ 
pect any neutrinos in a detector like Super-Kamiokande. 
Detections with Super-Kamiokande would be limited to 
the Milky Way and surrounding satellite galaxies. How¬ 
ever, with Hyper-Kamiokande, a future water-Cherenkov 
detector with a proposed volume ^20 times that of 
Super-Kamiokande, the expected number of neutrinos in 
the first ~10 ms would be 0(1) and maybe a few over 
the entire lifetime of the disk, thus detectable especially 
in light of further timing input provided by gravitational 


D. Magnetic Effects 

We have also investigated the possible impact of mag¬ 
netization within a timescale of ~ 15ms from the onset 
of the merger. To this end, we consider for simplicity 
the intermediate case (DD2) and endow the stars with 
a magnetic dipole of maximum strength 10 13 G. The 
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dynamics of the merger and early post-merger present 
several instabilities and processes capable of significantly 
amplifying the magnetization [20l l38l 97, 981 . Indeed, 
recent local simulations have shown that some of these 
instabilities can drive the magnetic energy density up to 
a level approaching equipartition with the kinetic energy 
of the turbulent fluid on a timescale given by the turn¬ 
over time of the resulting eddies |22l 199] . However, these 
process occur on very fine scales that require numerical 
resolution (and computational resources) generally be¬ 
yond the reach of codes such as ours. Furthermore, with 
grid-based codes fully capturing these dynamics require 
resolutions of the order of < 0.1m. 

Instead of such exorbitant resolutions, we introduce an 
“effective driver” mechanism (similar to those recently 
used in other works mm) to account for this am¬ 
plification by including a term in the evolution equation 
for the magnetic field of the form 

d t (V7 B i ) = ...+a^B i . (26) 

Here we have defined the amplification factor as 

£ = £o©(p[w* - uJt hresh \)g(t)B k v k , (27) 

with Q(x) the Heaviside function so that the amplifica¬ 
tion factor vanishes when the z-component of the fluid 
vorticity, uj z , is below some threshold <z^ hresh . This 
threshold is chosen to only activate this term for the ed¬ 
dies with wavelengths shorter than a third of the MNS 
radius. (This avoids, in particular, adding an additional 
increase to the magnetic field due to the star’s differen¬ 
tial rotation.) The factor B^v k restricts the amplifica¬ 
tion only to the components along the velocity lines that 
stretch the magnetic field. The function g(t) restricts 
when this sub-grid model is active as described below. 
We note that the amplification introduced by this term 
is not explicitly divergence free. However, we monitor 
the divergence of the magnetic field, and find that, while 
it increases at merger, the rate of increase relative to the 
rate of increase of the field strength itself is comparable 
to that seen when the sub-grid model is inactive. 

This added term must be adjusted to match the expec¬ 
tations established by local simulations [22] . The term 
should result in the amplification of the magnetic field 
during merger until the associated magnetic energy den¬ 
sity reaches strengths close to equipartition. In partic¬ 
ular, small-scale dynamos act to amplify the magnetic 
field until the magnetic energy density reaches roughly 
« 60% of the local turbulent kinetic energy density [22]. 
We therefore allow this sub-grid model to act for a few 
milliseconds until the maximum magnetic field reaches 
10 16 -10 17 G (see for instance Ref. |102| for a formal 
calculation of this turbulent kinetic energy). The bulk 
dynamics, combined with this sub-grid model, give rise 
to a MNS with a strong and largely toroidal magnetic 
field structure. It is worth mentioning that we have tried 
an amplification factor without the factor B^v k and ob¬ 
tained similar results. The ensuing topology is illustrated 
in Fig.[B| 


We display the magnetic energy of the star in the top 
frame of Fig. [l7| We call the evolution with no sub-grid 
model our low magnetization case and compare it to the 
case with the sub-grid model active, our high magneti¬ 
zation. Also included in Fig. [17] are plots of the maxi¬ 
mum density (middle panel) and the resulting GW signal 
(bottom panel). These figures reveal essentially no dif¬ 
ference between the two magnetization cases within this 
timescale, despite the high magnetization case reaching 
almost 10 17 G. Nevertheless, subtle differences arise with 
finer-scale phenomena. 

We display the density, electron fraction, and magnetic 
field strength for both the low and high magnetization 
cases on a meridional plane in Fig. [18] A difference in 
electron fraction appears in the polar regions, but other 
differences are quite difficult to distinguish. And so we 
compute the distributions of different quantities and dis¬ 
play them in histograms. In particular, Fig. [I9| displays 
the electron fraction for unbound material. While the 
low case has a flat distribution in Y e , the high magneti¬ 
zation case peaks at roughly Y e « 0.15. Note also that 
the high case ejects roughly twice the material of the low 
case, presumably due to the additional magnetic pres¬ 
sure. The angular distribution of the unbound matter is 
also shown in Fig. [19] The angular distribution suggests 
that the extra material unbound by the high magnetic 
field occurs primarily along equatorial directions. 

The velocity distributions for both bound and unbound 
matter are shown in Fig. [20] The bound matter distri¬ 
butions for the two cases are largely the same, but the 
unbound matter has a slightly higher average (but not 
peak) velocity for the high magnetization case. 

Magnetic fields can also redistribute angular momen¬ 
tum and transport it from the inner to outer regions of 
the MNS [ 103 . Generally the processes that do this (e.g. 
magnetic breaking) take place over a timescale given by 
t ~ Rmns / v a with va the Alfven speed and R the radius 
of the star. For fields amplified to strengths of « 10 16 G, 
the time scale is roughly r ~ 100 ms and so a complete 
exploration of this possibility is rather costly. Neverthe¬ 
less, we can extract indications of such a redistribution 
by comparing the behavior of the angular momentum 
for the low and high magnetization cases. Fig. 21 shows 
the radially integrated value of the specific, z-component, 
angular momentum, l z , within a given radius. Apparent 
in the figure is a small difference between the high and 
low cases. In particular, the highly magnetized case has 
less angular momentum at small radii than the weakly 
magnetized case, indicating that angular momentum has 
been transported outward. As a consequence, the con¬ 
tribution of rotational support against collapse weakens 
in the more magnetized case which could help induce a 
more prompt collapse to a black hole. 

We note that the atmosphere used in these two mag¬ 
netized cases is larger than the unmagnetized cases by 
two orders of magnitude. Because the differences that 
we see are largely where the outer material mixes with 
the atmosphere, it becomes difficult to isolate quantita- 
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tive physical effects from atmosphere effects. Thus, the 
effects just discussed should be taken as qualitative indi¬ 
cations while we work to further improve the treatment 
of the atmosphere. 



FIG. 16: Density (colormap) and magnetic field vectors (blue 
arrows) for the high magnetization case at t ~ 15ms after the 
merger, when the system has settled into a quasi-stationary 
configuration. The vectors show that the magnetic field is 
mostly toroidal. 


IV. CONCLUSIONS 

In this work we have explored the rich phenomenol¬ 
ogy of binary mergers of equal mass neutron stars using 
realistic equations of state, neutrino cooling, and elec¬ 
tromagnetic effects. Our results indicate a tight connec¬ 
tion between different possible observables and the neu¬ 
tron star EoS. In addition to the impact of the equation 
of state on gravitational waves that are emitted during 
and after the merger, both neutrino production and the 
properties of ejected material are strongly dependent on 
the EoS (which, in turn, would have a strong impact on 
electromagnetic signals from radioactive decay). We also 
developed a sub-grid model designed to capture the am¬ 
plification of the magnetic field generated by turbulent 
flow during the merger. We find indications of increased 
angular momentum transport with the amplified mag¬ 
netic field, though a detailed analysis is left to a future 
work. 

The EoS imprints subtle differences on the expected 
GW signals prior (but close) to merger, and these differ¬ 
ences become much more significant during merger and 
afterward. With softer equations of state, for example, 
the merger occurs at higher frequencies than for stiffer 
equations of state. The orbital frequency at merger ex¬ 
erts a strong influence on the resulting MNS. In partic¬ 
ular, an analysis of the dominant frequency components 
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FIG. 17: Comparison of low and high magnetized cases for 
the DD2 binary, (top) The total magnetic energy for the do¬ 
main. (middle) The maximum density for the two cases, (bot¬ 
tom) The primary mode of the GW signal. The merger takes 
place at 14.3 ms and during this merger the magnetic field 
grows significantly. Even for the high magnetization case in 
which B max ~ 10 17 G, the stellar dynamics remains largely un¬ 
affected as shown in the bottom two panels primarily because 
the magnetic field is mostly toroidal. Other quantities, such 
as the temperature and the neutrino radiation rate, also dis¬ 
play small differences between low and high magnetizations. 


of the gravitational signal after merger reveals that they 
can be roughly approximated by this orbital frequency 
at merger together with the characteristic oscillations of 
the MNS. Our results suggest that the three waveforms 
presented here may be distinguishable by aLIGO for op¬ 
timal configurations up to distances of « 120 Mpc (with¬ 
out assuming other techniques like stacking of signals to 
further enhance this possibility). These results are con¬ 
sistent with those found in [32] . 

An analysis of the neutrino production indicates that 
EoS yielding more compact stars (softer EoS) produce 
the largest neutrino luminosity with the highest aver¬ 
age neutrino energy. As discussed, this is intuitively ex¬ 
pected as the collision takes place deeper in the gravi¬ 
tational potential with a more violent collision. We es¬ 
timate that neutrinos produced by these mergers could 
be observed anywhere in the Milky Way or surrounding 
satellite galaxies with current neutrino detectors such as 
Super-Kamiokande. Using larger detector volumes like 
the proposed Hyper-Kamiokande detector, the mergers 
studied here predict that we would see 0(1) neutrinos 
within the first ^10-20 ms after merger. 

An analysis of the ejecta indicates that only the softest 
EoS yields sufficient material, and of the right quality, to 
power a kilonova event peaking in the infrared. Indica¬ 
tions such as this are of significant current interest in light 
of the recent observations reported in [I4j [15] and m- 
These observations found an infrared emission occurring 
roughly a week after sGRB 130603B and about thirteen 
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FIG. 18: Snapshots of density, electron fraction, and mag¬ 
netic field strength in the meridional plane x = 0 at t ~ 9ms 
after the merger, when the system has settled into a quasi¬ 
stationary configuration. The low magnetization case is dis¬ 
played on the left and the high magnetization case is shown 
on the right. 
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FIG. 20: Distributions of speeds for stellar material for the 
low/high magnetized cases at t ~ 9ms after the merger. 
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FIG. 21: Comparison of the radially integrated specific angu¬ 
lar momentum component l z for the low and high magnetized 
cases using the DD2 EOS at 3.2 and 8.2ms after the merger 
takes place (t — 17.5 and 22.5ms respectively). Angular mo¬ 
mentum transport in the highly magnetized case is stronger 
and, as a result, its integrated value in the central region of 
the MNS is smaller than its lower magnetized counterpart 
case. 



Y e polar angle [radians] 


FIG. 19: Distributions of Y e and angular direction for un¬ 
bound stellar material for the low/high magnetized cases at 
t ~ 9ms after the merger. 


days after GRB-060614, consistent with predictions of 
the decay of heavy r-process elements. That softer EoS 
merge through a more violent collision sets the proper 
stage for a large amount of quite neutron rich ejecta, 
a requirement for producing heavy elements through r- 
processes which, in turn, can produce an observable in¬ 
frared emission through their radioactive decay. 

Finally, it is then interesting to connect these obser¬ 
vations with our numerical results to speculate about 
the neutron star EoS. As has been discussed in, for ex¬ 
ample, Ref. [36], non-vacuum compact binaries are the 


prime candidates for producing heavy elements through 
r-process nucleosynthesis. Such binaries are favored be¬ 
cause of: (i) the mounting evidence connecting compact 
binary mergers with sGRBs (e.g. m) and (ii) the ob¬ 
servational indications [35; that supernovae fail to pro¬ 
duce such heavy elements. If indeed these binaries are 
the primary producers of the heavy elements and if 
sGRB 130603B and GRB 060614 are representative, then 
one can speculate that the neutron star EoS is a soft 
one. Simulations presented here and also in Ref. [33] in¬ 
dicate that intermediate and stiff EoS are unable to eject 
enough material with a sufficiently high neutron com¬ 
position, and this result, at the least, is suggestive that 
neutron stars are described by a soft EoS. 

However, this line of reasoning is speculative and fur¬ 
ther numerical work is needed. In particular, the proper¬ 
ties of ejecta resulting from the merger of unequal mass 
neutron stars are yet to be analyzed, but simulations 
have shown that mergers of unequal binaries can gener¬ 
ate accretion disks more massive than those produced by 
equal mass binaries nnanns]. It thus would appear pos¬ 
sible that the amount of ejecta and its properties might 
depend sensitively on the mass ratio. Such studies are 
thus an important priority, although observational evi- 
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dence points to small mass ratios in binary neutron star 
systems m- 

An alternative path could be provided by black hole- 
neutron star binaries. (Ejecta estimates in such merg¬ 
ers have been examined recently in ED HHH].) However 
for sufficient ejecta with the required neutron richness: 
the black hole spin has to be quite high, the ratio of 
black hole mass to neutron star mass should be small, 
and the stellar EoS needs to be stiff. These ingredients 
enhance the possibility of tidal disruption. However ob¬ 
served black hole masses and stellar black hole spin es¬ 
timates (see e.g. ESI), if they hold generally, indicate 
black hole masses > 8M© with any possible spin. There¬ 
fore, the likelihood of r-process nucleosynthesis occurring 
in black hole-neutron star binaries remains unclear. 

It is certainly quite interesting that these two pos¬ 
sible paths—binary neutron star or black-hole neutron 
star mergers—seem to prefer two very distinct possibil¬ 
ities for the EoS to yield kilonova events like the ones 
reported. Excitingly, based on current observational ev¬ 
idence, a single gravitational wave observation tied to a 
sGRB might be enough to single out the correct case. 


Appendix A: The primitive solver 

High-resolution shock-capturing schemes integrate the 
fluid equations in conservation form for the conserva¬ 
tive variables, while the fluid equations are written in 
a mixture of conserved and primitive variables. It is well 
known that the calculation of primitive variables from 
conserved variables for relativistic fluids requires solving 
a transcendental set of equations. Our method for solving 
these equations with a finite-temperature EoS is a mod¬ 
ification of the algorithm that we use for the ideal gas 
EoS; the most significant change being that the internal 
energy must be calculated separately from the pressure 
using the table. 

The evolved conserved variables are defined as 


D = pW (Al) 

Si = ( hW 2 + B 2 ) Vi - (B j Vj) Bi (A2) 

r = hW 2 + B 2 - P - I {[B^f + A.) (A3) 

DY e = pWY e . (A4) 

The dominant energy condition places constraints on the 
allowed values of the conserved variables 

D> 0, S 2 <{D + t ) 2 , DY e > 0, (A5) 

and depending on the EoS the second condition can be 
sharpened to S 2 < (2 D + r)r mo]. These constraints 
may be violated during the evolution due to numerical er¬ 
ror, and they are enforced before solving for the primitive 
variables. A minimum allowable value of the conserved 
density D v ac is chosen, and if D falls below this value 
we set v 1 = 0 and D —>• D vac . We choose D vac as low as 


possible for the unmagnetized neutron star binary, which 
is about nine orders of magnitude smaller than the ini¬ 
tial central density of the stars. For the magnetized case 
our solver fails for such tenuous atmospheres, and it is 
increased by two orders of magnitude. If the second in¬ 
equality is violated, then the magnitude of Si is rescaled 
to satisfy the inequality. Finally, DY e is required to sat¬ 
isfy the constraint on D, and the computed value of Y e 
must be in the equation of state table. 

The primitives to be found are the density, p, pressure 
P, electron fraction T e , internal energy density e (i.e., or 
the temperature T, once the EoS is known) and velocity 
three-vector v l . The magnetic field is at the same time 
conserved and primitive field. We write the transcenden¬ 
tal equations in terms of the new rescaled variable 


_ hW 2 
= 


(A6) 


where h is the total enthalpy h = p(l + e) + P, and 
calculate Y e from the evolution variables DY e /D. Fol¬ 
lowing m, we rescale the conserved fields in order to 
get order-unity quantities, namely 

q = t/D, r = S 2 /D 2 , s = B 2 /D, t = BiS \/ D 3/2 . 

(A7) 

Then, using data from the previous time step to calculate 
an initial guess for x, we iteratively solve these equations 
for x within the bounds 


1 -\- q — s>x>2-\-2q — s , (A8) 


so that the final procedure can be written as 

1. From the equation for S l Si , calculate an approxi¬ 
mate Lorentz factor W 


x 2 r + (2x + s) t 2 
x 2 (x + s) 2 


2. From the definition of D, calculate 


P = 


D 

W 


3. From the definition of r and the total enthalpy, 
calculate 


emW-1 


+ 4- (l-W 2 ) 
VE V > 


2x 2 + 2W 2 


4. Use the EoS table to calculate the pressure 

P(p,e,Y e ). 

5. Update the guess for x by solving the equation 
f{pc) = 0 using the Brent method, being f{x) just 
the definition of the unknown x, 


f(x) = x- 


l + e + 


P(p,e,Y e ) 


W 
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The root of f(x) = 0 from Step 5 becomes the new guess 
for x, and this process is repeated iteratively until the 
solution for x converges to a specified tolerance, which is 
ensured if there is a physical solution within the bounds. 
One advantage of this algorithm is that f{x) is a function 
of a single variable, and, in contrast to a multiple variable 
search for a root, robust methods can be used to find any 
root that can be bracketed. 

Because of numerical error, a solution to these equa¬ 
tions may either fall outside the physical range for the 
primitive variables, or a real solution for x may not ex¬ 
ist. The solutions for p, T, and Y e are, at a minimum, 
restricted to values in the table, and they are reset to new 
values (the minimum allowed value plus ten percent) if 
necessary. If a real solution for the primitive variables 
does not exist, the primitive variables are interpolated 
from neighboring points, and the conserved variables are 
reset to be consistent. If a valid interpolation stencil can 
not be constructed because the solver also failed at the 
neighboring points, then the update fails, and the run is 
terminated. This failure occurs very rarely and may be 
remedied by slightly increasing the density floor D vac . 

Appendix B: Convergence tests 

The inspiral and merger of two neutron stars is a signif¬ 
icant test due to the inherent asymmetry of the problem 
that tends to “excite” many, if not all, terms in the equa¬ 
tions. Resolving both the motion of two compact objects 
as well as the large gradients at their surfaces require sig¬ 
nificant resources. The merger itself is a very dynamic 
process with a large range of densities. Therefore, a study 
of the convergence of the numerical solution is necessary 
to assess the validity of our results. 

We have evolved a binary using the DD EoS with three 
different resolutions (i.e., low, medium, and high) corre¬ 
sponding to Ax = {275, 230, 192} m. The stellar sep¬ 
aration and the total baryonic mass, which should be 
strictly conserved during the evolution, are displayed in 
Fig. [22] Interestingly, the merger happens at earlier times 
for higher resolutions, a trend opposite of that reported 
elsewhere. This might indicate that the source of our 
numerical errors has a different sign. In any case, the 
results converge at order ^1.9. 

The convergence of the principal l = m m 2 mode 
of T 4 is shown in Fig. [23] Finally, the distribution of 
the velocity of the fluid, computed just few milliseconds 
after merger, is displayed for the different resolutions in 
Fig.m 
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were performed at XSEDE and Scinet. 


FIG. 24: Histogram of the velocities, computed a few mil¬ 
liseconds after merger for the same cases as shown in Fig. |22| 
Notice that the medium and high resolution results have 
largely converged. 
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